Fluctuation-induced collective motion: A single-particle density analysis 
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In a system of noisy self-propelled particles with interactions that favor directional alignment, 
collective motion will appear if the density of particles increases beyond a certain threshold. In this 
paper, we argue that such a threshold may depend also on the profiles of the perturbation in the 
particle directions. Specifically, we perform mean-field, linear stability, perturbative and numerical 
analyses on an approximated form of the Fokker-Planck equation describing the system. We find 
that if an angular perturbation to an initially homogeneous system is large in magnitude and highly 
localized in space, it will be amplified and thus serves as an indication of the onset of collective 
motion. Our results also demonstrate that high particle speed promotes collective motion. 
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The interesting phenomena of flocking in animals [1| 
and self-organized patterns in motile cells 041] are cur- 
rently driving the intense theoretical study of collective 
motion among self-propelled particles |7H17| . In particu- 
lar, a comprehensive linear stability analysis on the onset 
of collective motion from the perspective of Boltzmann 
equation has recently appeared [lg]. Models for collec- 
tive motion usually involve motile particles that possess 
alignment interactions and angular noise. Collective mo- 
tion is then observed if the density of particles increases 
beyond a certain threshold. Besides density fluctuations, 
fluctuations in the heading directions of the particles con- 
stitute another important aspect of the system. Here, we 
study a minimal model for collective motion and show 
that the threshold for collective motion transition may 
depend on the profiles of directional fluctuations. Specif- 
ically, we find that an initial directional perturbation to 
a spatially homogeneous system will be amplified, if the 
perturbation is large in magnitude and is highly localized 
in space. We also demonstrate that high particle speed 
promotes collective motion. 



To achieve our results, we first write down the Fokker- 
Planck equation describing the single-particle density dis- 
tribution of the system in Sect. [II] We then investigate in 
Sect. IIIII the equation in the Fourier space of the direc- 
tional component, and argue that only the lower order 
modes are important at the onset of collective motion. 
As a result, the dynamics of the distribution function can 
be captured by a set of three nonlinear coupled differen- 
tial equations, which we subsequently study with linear 
stability analysis in Sect. IIVI In Sect. [V] we go beyond 
the linear stability regime by investigating the dynamical 
equations perturbatively and numerically. 



In this work, we follow [16| and consider a minimal 
model for collective motion in two dimensions, where ev- 
ery particle is assumed to have constant speed and that 
their interactions consist only of directional alignment 
mechanism. Noise is incorporated in the direction of 
travel. Specifically, let there be N particles in a volume 
of V, their equations of motion are: 



dr; 2u 

= — v(0i) 
at 7t 



(1) 

/ 2Dr h (t) (2) 

where 1 < i < N, R = (n, . . . , r N ), 9 = [9 U 6 N ), 
v(0) = (cos 0, sin 0), and the noise is assumed to be Gaus- 
sian characterized by the following moments: 

(m(t)) = o , ( ni (t) Vj (t')) = s tJ s(t-t') (3) 

Moreover, the alignment interaction is assumed to be of 
very short range and thus can be approximated by a delta 
function: 



U(R,e) = --Y <5 (2) (r 4 -r,)cos( 



■ (4) 



i<j 



If we denote the probability distribution of the density of 
particles in the state (R, 6) at time t by f(t, R, 0), then 
the Fokker-Planck equation corresponding to the system 
is El: 



df 
dt 



9 9 

i<j 



d (2 \r l ~r 3 ) S m(d. l -d 3 )f] . (5) 



Naturally, we are not interested in all the information 
captured by /, and we will from now on focus on the 
single-particle density function, p, 
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Prom Eq. ((5|), we can express p in terms of the two- 
particle density function p^: 



dp(r,8) =D d 2 p(r,6) 



dt 
where 



dO 2 



2 u 



cos e d -£m +si n« dpir > 9) 



dx 



dy 



p^{r 1 ,9 1 ,r 2 ,9 2 



do' I &r'8 {2) {Y-Y')sin{e-e') P {2) (Y,e,r',e') 

_ (M) J dr 3 • • • dr N d0 3 ■ ■ ■ d9 N f(R, 9) 



(7) 



(iV-2)! 



(8) 

The above manipulation is akin to the BBGKY hier- 
archy formalism [2(ij |. To continue with our analyti- 
cal treatment, we will ignore the second ordered cor- 
relation and adopt the product distribution assumption: 
p^(r, 9, r', 9') = p(r, 9)p(r', 9'). This assumption is sim- 
ilar to the molecular chaos assumption in the context of 
Boltzmann equation, and is also adopted in [TEl fTsj j . 

By Fourier transforming the above equation with re- 
spect to the angular variable, 9, we have: 

d t pn(r) = -Dn 2 p n (r) - u d x (p n+1 {r) + p n -i(r)) 



+idy(p n - 1 (r) - /5„+i(r)) 



-gn 



/3_i(r)p n+ i(r) - pi(r)p„_i(r) (9) 



where p(r, 9) = EZ-oo Pn(r)e~ ine . 

Since /5„(r) are complex, we will denote them by 
o„(r) + ib n (r) where a n and b n are real functions. In 
relation to the original density function, we have 



p(r,6) = ^K(r)+i6 n (r)]e- ine (10) 
r) + 2 ^ ^ri(r) cos(n9) + b n (r) sin(n9) 



ao(r) 



n>l 



where for the second equality, the following conditions 
for the a n and b n have been employed: 



(11) 



which are due the fact that p is real. Writing Eq. |H] in 
terms of the a n and b n , we have for n 6 Z, 



dta n — —Dn 2 a n — u 
d t b n = ~Dn 2 b n - u 



d x (a n+ i + a„_i) - dyibn-x - b n+1 ) - gn ai(a n+i - a„_ x ) + 6i(6„_i + b n+1 ) (12) 
d x (b n+1 + b n -i) + d y (a n+1 - a n _i) - gn ai(6„+i - 6„_i) - &i(a„+i + a„_i) . (13) 
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Note that the arguments (i, r) in a„ and b n are omitted 
in the above equations to ease notation. 



III. MEAN-FIELD APPROXIMATION 

To avoid having to deal with the above infinite set of 
differential equations, we will sort to truncate the number 
of differential equations to be considered. To do so, we 
first study the system in a mean-field manner plj . i.e., 
we set p n (t,r) — p n (t) for all r. Eqs (TT^j) and (|13l) then 
become 



da n 
dt 

db n 
dt 



-Dn 2 a„ 



gn 



ai(a n 



+i 



+6i(6„_i + b n+ i) 



-Dn 2 b rj 



gn 



ai(b 



n+l 



-bi(a n+ i + a„_i) 



(14) 



(15) 



Note that dao/dt = due to the fact that a corre- 
sponds to the overall density of the system, which does 
not change. 



Let us assume that the b modes are not excited at t = 
and so we need only focus on the a modes. By inspecting 
Eq. [10l we see that the omission of the b modes is the 
same as focusing only on angular perturbation of the form 
cos(#), i.e., the particles are more likely to be heading in 
the positive x direction. With this simplification, the 
first three modes are of the form: 



dai 
~dT 
da 2 
~dT 
da 3 
dt 



gaoai — (Dai + gaia 2 ) (16) 
2ga\ - (4Da 2 + ga ia3 ) (17) 
iga\a 2 — {QDa 3 + ga\a^) . (18) 



At the onset of collective motion (CM) from a spa- 
tially and angularly homogeneous system, we expect that 
\a n \ <ti 1 for n > 1. Let us define e as max„>i \a n \ at the 
onset of CM, we see that only da\/dt is of order e while 
all the time-derivatives for the higher order modes are of 
order e 2 . Furthermore, the coefficients associated with 
the damping term D for the n— th modes scale with n 2 , 
which further suggests that only the lower order modes 
are imp ortant. Another corroborating evidence is from 
[l5l Il8| where the authors employed their scaling ansatz, 
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which is supported by their numerical simulations, to ar- 
gue that the higher order modes are indeed negligible 
at the onset of CM. Based on all these reasons, we will 
truncate the original dynamical equations, Eq. ([§]), by 
omitting all p n for n > 2. Focusing again only on the a 
modes, we have 



dap 
dt 
dai 

~d7 

da 2 



= 



(-D + ga Q - ga 2 )a 1 
—ADa 2 + 2ga\ . 



(19) 
(20) 
(21) 



A simple fixed point analysis on the above equations in- 
dicate that the existence of non-zero fixed point for a\ 
and ci2 is only possible when 



ga -D>0 



(22) 



This condition has previously been derived in jig ]. Ex- 
pectedly, the above condition indicates that collective 
motion is facilitated by having strong interaction (g), 
high particle density (ao) and weak noise (D). 



FIG. 1: The density profiles of q, /3 and 7 at different 
times obtained by numerically solving the set of differential 
equations in Eqs (|30|) to (|32[) . with the following parameters: 
u = 1/a/3, and D = g = (, = a = 0.1. 



1.5 



0.5 



0.4 



0.2 



-0.2 









✓ " 
■ ■ / ■' 
















/ 












-0.5 









0.5 




















t = 
















t = 0.1 






// 










t = 0.35 






\* i 

rJ 

- '_L 




\* 

\y 
\\ 

. ■ v^s 










-0.5 






/ 

/ 




0.5 












/ 












-0.5 






X 




0.5 







IV. LINEAR STABILITY ANALYSIS 

We now continue with our truncation approximation, 
but with the spatial variable re-installed into Eqs (fT4|) 
and (TJ5J) . Before we start to analyze the set of differen- 
tial equations, we note that by inspecting Eqs (I12I) and 
(p~3|) . we see that the a n >i and 6„>i modes are coupled 
exclusively to different spatial dimensions - the x and y 
dimensions respectively. In other words, if the system 
is initially homogeneous in the x dimension, then the x 
dimension will remain homogeneous, and vice versa. We 
will therefore, as in the previous section, assume that the 
b modes are not excited and focus only on the a modes. 
With this simplification, we arrive at the following dy- 
namical equations: 



dt.a = —2ud x /3 

d t P = -D/3~ud x (a 



- 2gP 2 , 



7) 



(23) 
(24) 
(25) 



where we have used the Greek letters a, /3 and 7 to denote 
ao, ai and 0,2 respectively. 

The fixed point in the homogeneous phase corresponds 
to a = 1, /3 = and 7 = where we have set the unit 
length in such a way that the density of the particles is 
one. We now perform linear stability analysis on this 
fixed-point by considering the linear response of the sys- 
tem to a small perturbation of the form: 



a 

7 



1 + Ac xt+iqv 
Be xt+iqy 



(26) 
(27) 
(28) 



where A,B,C <C 1 and q is an arbitrary frequency. Sub- 
stituting the above into Eqs (|23p and gives the fol- 
lowing condition on A: 



A(A + D - g){\ + 4D) 
3A + 8D 



2 2 
-u q 



(29) 



which indicates that A > if and only if g > D. In 
other words, we have recovered the condition found in 
our previous mean-field analysis (c.f. Eq. (|22"]) ). 

Although this result is consistent with what we found 
in the previous section, pieces of the picture at the onset 
of CM are still lacking. For instance, the phase transition 
condition found here does not depend on the speed of the 
particle u. This is unsatisfactory because we know that 
long range order would not be possible if it = [a Q ■ 
Moreover, it is desirable to see how the coupling between 
the spatial and temporal dimension affects the rise of 
the excited mode (3. To gain insight in these questions, 
we will go beyond the linear stability regime and analyse 
the dynamical equations with perturbative method in the 
next section. 



BEYOND LINEAR STABILITY ANALYSIS 



In this section, we will study Eqs (|23|) to (f25|) perturba- 
tively. Specifically, we will assume that D,g <C 1, in the 
units of distance and time set by having a(t = 0, x) = 1 
and u = 1/V3. Physically, these assumptions amount to 
limiting our discussion to the regime where the particles' 
angular fluctuations and interaction strength are small. 
Furthermore, since we are primarily interested in the dy- 
namics at the onset of CM, we will assume that g/D is 



4 



of order unity (c.f. Eq. (1221) ). These assumptions allow 
us to employ D (or equivalently, g) as the expansion pa- 
rameter in our perturbative treatment. In contrast to the 
linear stability analysis in the previous section where the 
perturbation magnitude is assumed to be small enough 
that the nonlinear term is negligible, the perturbative ap- 
proach adopted here allows us to study the effects of the 
nonlinear term on the dynamics. 

In the aforementioned units, Eqs (|23|) to (1251 are: 



d t a 



V3 



(30) 



d t p = -Dfi - ~jf x {a + 7 ) + gp(a - 7 ) (31) 



d a = -4^7 - ±=d w p + 2gp 2 . 

We now expand a as: 

a = a + Dm + 0(D 2 ) , 



(32) 



(33) 



and similarly for /3 and 7 . The zero-th order (in D) terms 
follow the following differential equations: 



FIG. 2: The temporal evolutions of (a) A (c.f. Eq. ((45j)) 
and (b) B (c.f. Eq. (JUJ), obtained by numerically solving 
the set of differential equations in Eqs (|30[) and (|32|l . with 
D — g = £ = 0.1. (c) The zoom-in plot of B(t) at small 
time with the three curves corresponding to the theoretical 
expressions given in Eq. (149 1 . 
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d t a 



V3 



d x /3o 



dtPo = --j=d x (ao+'yo) 
dtlo = --^=d x j3o . 



(34) 
(35) 
(36) 



The above set of differential equations can be solved by 
employing the Laplace-Fourier Transform method. For 
the initial conditions of ao(t = 0,x) = 1, /3o(t = 0,x) = 
£ e ~ y 2 /(2a 2 ) / y /2^ a and 7o ( f _ QjX ) = o. The solutions 
are: 



«o 



1 



wherc 



Po = 



7o = 



2V27r(T 



6it<j 

U~ +U + 



2V67TCT 



exp 



(x ± tf 
2a 2 



(37) 
(38) 
(39) 

(40) 



In other words, a Gaussian perturbation in /3 at t = 
splits into two Gaussian distributions traveling in oppo- 
site directions with unit speed (as a result of setting u 
to I/a/3). The perturbation also induces in a and 7 two 
solitary waves in the form a Gaussian distribution trav- 
eling with unit speed in the positive direction, and an 
inverted Gaussian density wave traveling in the opposite 
direction (c.f. Fig. [T]). This is akin to the stripe traveling 
wave pattern found in the CM phase [1(1 [3 ■ 



The differential equations governing the first-order 
terms are: 



<9 t ai = --j=d x pi 



(41) 



dtPx = -/?o--^^(a 1 +7i) + ^/3o(ao-7o)(42) 



(43) 



dsn = -4 7o - -j=d*lh + 



where ao, A) and 70 above are now given by Eqs (|37|) to 
(|39p . The initial conditions for the above equations are: 
ai(t = 0,x) = Pi(t = 0,x) = 71 (t = 0,x) = 0. 

Before we attempt to study the above set of differen- 
tial equations, let us look for some meaningful quantities 
that quantifies the effect of the initial perturbation. For 
instance, the temporal evolution of total increase in /3\ 
due to the initial perturbation can be obtained from Eq. 
63): 

/ r 00 \ r°° 



(44) 



This indicates that when summed over the whole space, 
the P mode is amplified only if g > D. This is again 
consistent with the result we obtained in Sect. IIIII and 
Sect. [TVl 

Besides the above quantity, the following two quanti- 
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FIG. 3: The temporal evolutions of (a) A (c.f. Eq. J45])) 
and (b) B (c.f. Eq. (|45jl). obtained by numerically solving 
the set of differential equations in Eqs (|30[) an d (|32|l . with 
g = b — a = 0.1. The zoom-in plot of Bit) at small time with 
the three curves corresponding to the theoretical expressions 
given in Eq. J49j. 



FIG. 4: The evolution of B(t) obtained by numerical sim- 
ulating the differential equations Eqs (|12p and (|13[l with the 
initial condition discussed in Sect. [V] The parameters are 
g = D = £ = a = 0.1. In the simulations, only the 0-th to the 
m-th a modes are included. In other words, we assume that 
an>m = and fe„gz = 0. 
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A(t) 



dxa(t, x) ; B{t) 



dx[/3(t,x)-0(O,x)] . 

(45) 

Namely, A and B correspond to the responses in the den- 
sity (a) and in the f3 mode in the direction of the angular 
perturbation. These are in fact arguably better quanti- 
ties to consider as they capture the directional nature of 
the perturbation. From Eq. (I42I) . we have: 



dB(t) _ Z(g-D) 



9^ 



At 



erffl 
37TCT \<y 

+-=[a>x(t,x = 0) +ji{t,x = 0)] 



V3 

In the appendix, we demonstrate that 



a 1 (t,y = 0)+j 1 (t,y = 0) 



9 e 



t + 0(t 3 



In other words, up to order 0(t 3 ), we have 



B(t) 



b(g-D)t , 5gt 2 



8\/37r \ o~ 



(46) 
(47) 

(48) 

(49) 



The third term in the right hand side above highlights 
the importance of the term £/<r, especially when D ~ g. 
Fig. [SJa) and (b) display the temporal evolutions of A(t) 
and B(t) by solving Eqs (f30|) and (|3~T1) numerically in the 
case of D = g. They clearly show the amplification of 
the initial perturbation, which we have taken as an indi- 
cation for the onset of collective motion in longer time. 
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Fig. Oc) demonstrates that at short time, the dynam- 
ics is well described by the expression in Eq. (|49p. Fur- 
thermore, due to the positive second term in the R.H.S. 
above, the formula for B{t) suggests that there is a pos- 
sibility of perturbation amplification even if D > g, e.g., 
when £/cr ^> 1. This is indeed shown to be the case in 
Fig. (3{a) and (b) , where amplification of the perturbation 
is seen for D/g ~ 1.1. In other words, a sharp pertur- 
bation in the angular domain is able to induce collective 
motion even if the density is below the phase transition 
threshold as obtained in the mean-field model. 

If we now restore the speed, u, and the initial density, 
c, where c = a(t = 0, x), into Eq. (|49|) . we have 



Bit) 



bc(cg — D)t bguct 2 ( £ 
2 + 8tt 



(50) 



Note that the speed of the particles only appears in the 
second term above, which is positive. Hence, the above 
formula suggests that the particle speed has a net effect of 
amplifying the initial perturbation and thus facilitating 
collective motion transition. 

As a verification on the validity of our approximation 
adopted in our analytical calculations, we numerically 
simulate Eqs (fT2|) and (fl3)) with higher order modes in- 
cluded and find that there are no discernible differences 
for the parameter range investigated in this work (c.f. 
Fig. S]). 



VI. CONCLUSION 

In summary, starting with a Fokkcr-Planck equation 
for a minimal model of CM, we derived a set of three 
coupled differential equations that describes the system 
at the onset of CM. We then studied the equations with 
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mean-field, linear stability, perturbative and numerical 
analyses, and found that if an angular perturbation is 
large in magnitude and highly localized in space, it will 
be amplified and thus serves as an indication of the onset 
of collective motion. Our calculations also demonstrate 
the importance of particle speed for collection motion 
transition. As a result, it is indicative that the critical 
point for CM may depend on the speed u, the pertur- 
bation magnitude b and the perturbation wavelength a. 
This is in contrast to the mean-field and linear stabil- 
ity analyses where only the hydrodynamic, or infinite- 
wavelength, mode dictates the onset of CM. Our results 
therefore highlights the importance of incorporating the 
nonlinear term into the analysis. 

The main limitation of this work is on the approxi- 
mation adopted - the omissions of higher order modes. 
While we believe that such an approximation is appro- 
priate at the onset of CM, it would be highly desirable to 
have a systematic method to incorporate the higher or- 
der modes into the dynamics. Besides the consideration 
of the higher order modes, singular perturbation method 
would also be needed to investigate the long-time be- 



haviour of the system [22| . We believe that these aspects 
would constitute two promising directions for future in- 
vestigation. 



Appendix 

We are unable to solve the set of differential equations 
shown in Eqs (j^Tj) to P5)l analytically. But since only the 
leading orders in x and t are of interests, we will replace 
the C/ ± in a ,A),7o (cf. Eqs (|3TJ) to ((55)1) by 



(x±tf 
2a 2 



{x±tf 



2a 2 



(A.l) 



With this simplification, the differential equations can be 
solved by the Laplace- Fourier Transform method and the 
relevant results are: 

ai(i,y = 0) =0(t 3 ) , <yi(t,y = 0) = ^t + 0(t 3 ) . 

7T(7 

(A.2) 
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